package twoD.IsoPfem.Plasticity.J2PlaneStress;

import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import oneD.fem.algorithm.ITensor;
import oneD.fem.algorithm.Tensor2;
import oneD.fem.algorithm.Tensor4;
import twoD.IsoPfem.CST;
import twoD.IsoPfem.IGeometry;
import twoD.IsoPfem.RealTable;
import twoD.IsoPfem.Plasticity.IHardeningLaw;
import twoD.IsoPfem.Plasticity.IPlasticElement;

public class PlasticCST extends CST implements IPlasticElement {

	private static final double TOL = 1e-6;

	private Tensor4 tangentOperator_;
	private Tensor2 sigma_;
	private ArrayList<PlasticState> state_;

	public PlasticCST(RealTable realTable, IGeometry T) {
		super(realTable, T);
		state_ = new ArrayList<PlasticState>();
		tangentOperator_ = new Tensor4(2);
		sigma_ = new Tensor2(2);
	}

	@Override
	public IVector computeSigma(double... xi) {
		return sigma_.getVectorForm();
	}

	@Override
	public IMatrix getCMatrix() {
		IMatrix C = new Array2DMatrix(3, 3);
		C.set(0, 0, tangentOperator_.get(0, 0, 0, 0));
		C.set(0, 1, tangentOperator_.get(0, 0, 1, 1));
		C.set(0, 2, tangentOperator_.get(0, 0, 0, 1));
		C.set(1, 0, tangentOperator_.get(1, 1, 0, 0));
		C.set(1, 1, tangentOperator_.get(1, 1, 1, 1));
		C.set(1, 2, tangentOperator_.get(1, 1, 0, 1));
		C.set(2, 0, tangentOperator_.get(0, 1, 0, 0));
		C.set(2, 1, tangentOperator_.get(0, 1, 1, 1));
		C.set(2, 2, tangentOperator_.get(0, 1, 0, 1));
		return C;

	}

	public void updateInternalState(vonMisesYieldCriterion yieldCriterion) {
		// trial state
		PlasticState current_state = state_.get(getCurrentStep());
		Tensor2 epsilon = new Tensor2(computeEpsilon());
		Tensor2 e = epsilon.getDeviatoricTensor();
		Tensor2 ep = current_state.getEpsilon_p().getDeviatoricTensor();
		Tensor2 s = e.add(-1, ep);
		final double nu = getRealTable().getNu();
		final double K = getRealTable().getK();
		BLAM.scale(2 * nu, s.getA());
		Tensor2 xi = s.add(-1, current_state.getBeta());
		double norm_xi = xi.dot(xi);

		double f = yieldCriterion.f(xi, current_state.getAlpha());
		if (f <= 0) {
			PlasticState p = new PlasticState(epsilon,
					current_state.getEpsilon_p(), current_state.getAlpha(),
					current_state.getBeta());
			state_.add(p);
			sigma_ = s.add(K, new Tensor2(s.dimension()));
		} else {
			// return-map
			final double sqrt_2_3 = Math.sqrt(2.0 / 3);
			double delta_lambda = computeDeltaLambda(xi, yieldCriterion);
			Tensor2 N = (Tensor2) xi.clone();
			BLAM.scale(1.0 / Math.sqrt(norm_xi), N.getA());
			double alpha = current_state.getAlpha() + sqrt_2_3 * delta_lambda;
			double tmp = sqrt_2_3
					* (yieldCriterion.getKinematicHardeningLaw().valueAt(alpha) - yieldCriterion
							.getKinematicHardeningLaw().valueAt(
									current_state.getAlpha()));
			Tensor2 beta = current_state.getBeta().add(tmp, N);
			Tensor2 epsilon_p = current_state.getEpsilon_p().add(delta_lambda,
					N);
			sigma_ = s.add(K, new Tensor2(s.dimension())).add(
					-2 * nu * delta_lambda, N);
			double theta = 1 - 2 * nu * delta_lambda / norm_xi;
			double theta_bar = 1
					/ (1 + (yieldCriterion.getIsotropicHardeningLaw().valueAt(
							alpha) + yieldCriterion.getKinematicHardeningLaw()
							.valueAt(alpha)) / (3 * nu)) - (1 - theta);
			Tensor2 I = new Tensor2(s.dimension());
			Tensor4 C = new Tensor4(s.dimension());
			C.scale(2 * nu * theta);
			C.add(K - 2 * nu * theta / 3, I.product(I));
			C.add(-2 * nu * theta_bar, N.product(N));
			tangentOperator_ = C;
			PlasticState p = new PlasticState(epsilon, epsilon_p, alpha, beta);
			state_.add(p);
		}
	}

	public int getCurrentStep() {
		return state_.size() - 1;
	}

	private double computeDeltaLambda(ITensor xi,
			vonMisesYieldCriterion yieldCriterion) {
		double deltaLambda = 0;
		double alphan = state_.get(getCurrentStep()).getAlpha();
		double alpha = alphan;
		IHardeningLaw H = yieldCriterion.getKinematicHardeningLaw();
		IHardeningLaw K = yieldCriterion.getIsotropicHardeningLaw();
		final double sqrt_2_3 = Math.sqrt(2.0 / 3);
		double nu = getRealTable().getNu();
		double eps = 1.0;
		while (eps > TOL) {
			double dH = H.valueAt(alpha) - H.valueAt(alphan);
			double g = Math.sqrt(xi.dot(xi)) - sqrt_2_3 * K.valueAt(alpha)
					- (2 * nu * deltaLambda + sqrt_2_3 * dH);
			double dg = -2
					* nu
					* (1 + (H.derivativeAt(alpha) + K.derivativeAt(alpha)
							/ (3 * nu)));
			deltaLambda -= g / dg;
			alpha += sqrt_2_3 * deltaLambda;
			eps = g;
		}
		return deltaLambda;
	}
}
